Effective simulations of interacting active droplets

Droplets form a cornerstone of the spatiotemporal organization of biomolecules in cells. These droplets are controlled using physical processes like chemical reactions and imposed gradients, which are costly to simulate using traditional approaches, like solving the Cahn–Hilliard equation. To overcome this challenge, we here present an alternative, efficient method. The main idea is to focus on the relevant degrees of freedom, like droplet positions and sizes. We derive dynamical equations for these quantities using approximate analytical solutions obtained from a sharp interface limit and linearized equations in the bulk phases. We verify our method against fully-resolved simulations and show that it can describe interacting droplets under the influence of chemical reactions and external gradients using only a fraction of the computational costs of traditional methods. Our method can be extended to include other processes in the future and will thus serve as a relevant platform for understanding the dynamics of droplets in cells.

www.nature.com/scientificreports/ Continuous theory of phase separation. We consider an isothermal, incompressible fluid in a closed system of volume V tot that consists of solvent and droplet material. The composition is described by the volume fraction φ(x, t) of the droplet material, while the solvent volume fraction is given by 1 − φ . The thermodynamic state of the system is governed by the free energy functional F[φ] = f (φ) + κ 2 |∇φ| 2 dV ; see ref. 28 , where κ is a parameter related to the interfacial tension, while phase separation is promoted by the free energy density f (φ) . For simplicity, we here consider a polynomial form, where the minima φ 0 in and φ 0 out are the equilibrium concentrations in a thermodynamically large system and b denotes the energy scale. The dynamics of the system follow from the continuity equation, ∂ t φ + ∇ · j = s . Here, j denotes spatial fluxes, which for simplicity will only be driven by diffusive processes, so hydrodynamic fluxes are neglected. Conversely, s is a source term related to chemical reactions that convert droplet material into solvent and vice versa 39 . Chemical reactions are typically local and often described by rate laws that depend on composition. Conversely, the non-local diffusive fluxes j are driven by gradients in chemical potential where ν is the molecular volume of the droplet material. Linear non-equilibrium thermodynamics implies j = −�(φ)∇µ , where �(φ) is a positive mobility 40 . Hence, is a fourth-order, non-linear partial differential equation requiring two boundary conditions. We here focus on the typical choice of no-flux conditions ( n · ∇µ = 0 ) and that solvent and droplet material interact identically with the system's boundaries ( n · ∇φ = 0 ), where n denotes the normal vector at the boundary.
In the case without chemical reactions ( s = 0 ), Eq. (3) reduces to the seminal Cahn-Hilliard equation 28 , which describes passive phase separation. In particular, two bulk phases with composition φ 0 in and φ 0 out typically emerge. These phases are separated by an interface of width w = 2 √ κ/b with surface tension γ = 1 6 √ bκ 27,39 . When chemical reactions are weak, this general structure is typically preserved, although long-term dynamics, like Ostwald ripening, can be strongly modified 39 . Strong chemical reactions can actually destroy droplets 39 and they might also lead to more complicated patterns 41,42 , which go beyond the scope of this paper. Instead, we here focus on situations where well-defined droplets with a thin interface are typical.
The dynamics given by Eq. (3) adequately describe phase separation, but it can be prohibitively costly to simulate due to multiple reasons: (1) The interface needs to be resolved, implying discretizations on the order of the typically small interface width w. (2) The equation contains fourth-order derivatives in space, which often limits the time steps. (3) Interesting dynamics often take place on very long time scales. For instance, during Ostwald ripening 43,44 , length scales in the system evolve as t 1/3 , requiring long simulations to capture relevant behaviour. However, since we are interested primarily in modelling dynamics of droplets, we circumvent these problems by focusing only on phase separation inside the nucleation and growth regime of the free energy density, and assume that sufficiently finite perturbations have already nucleated droplets. We employ the thininterface approximation 39 , which is a coarse-grained analytical formulation of the continuous model Eq. (3) valid when the system is subject to strong phase separation, low variation of volume fractions in the droplet phase and the dilute phase, and large droplet sizes compared to the interface width. This analytical approach was utilized earlier to study kinetics of many-droplet systems and effects of chemical reactions on such systems 39 .
In the next section, we elaborate on using the thin-interface approximation to build an effective droplet model describing the dynamics of droplets and dynamics of the dilute phase separately, instead of the full volume fraction field from the continuous model, thus effectively 'de-coupling' the description of phase separated droplets from the dilute phase.

Effective description of isolated droplets.
To build the effective model, we next derive approximate descriptions of the dynamics of the radius R and position x of an isolated droplet. Since we only consider spherical droplets with a thin interface ( R ≫ w ) and weak chemical reactions, we can use basic thermodynamics to derive the equilibrium concentrations inside and outside of the interface of the droplet, denoted by φ eq in and φ eq out , respectively. Due to surface tension effects, they are slightly elevated above the basal values φ 0 in and φ 0 out prescribed by the free energy density; see Eq. (1). To first order in the curvature of the surface, we have where l γ ,out and l γ ,in are capillary lengths 27 . In our case they read l γ ,in = ( and in three dimensions. The dynamics of the volume fraction field φ in inside the droplet is in principle described by Eq. (3), but since the composition typically hardly varies, we can linearize φ in around φ 0 in to obtain www.nature.com/scientificreports/ where D in = �(φ 0 in ) b is the diffusivity and k in = −s ′ (φ 0 in ) denotes the reaction rate 27 . Generally, positive rates ( k in > 0 ) stabilize the volume fraction φ in , while negative rates might destabilize it. However, the instability is suppressed when the droplet radius R is small compared to the reaction-diffusion length scale, ξ in = √ D in /|k in | 27 . Since we here consider weak chemical reactions, ξ in will be large, and we thus assume R ≪ ξ in in the following. We use this to solve Eq. (5) in stationary state in a system with angular symmetry using the boundary conditions φ in (R) = φ eq in and ∂ r φ in (0) = 0 . The analytical result allows us to estimate the diffusive flux j in inside the interface, where d is the space dimension; see Supporting Information, Section III. Production of droplet material inside the droplet ( s(φ eq in ) > 0 ) leads to an outward flux j in · n > 0 , which can drive droplet growth. Conversely, destroying droplet material ( s(φ eq in ) < 0 ) promotes shrinking droplets. Droplets might also grow if they take up material from the surrounding. Similarly to inside of droplets, the volume fraction φ out will typically vary only little, so we can linearize around the base value φ 0 out to obtain the reaction-diffusion equation where D out = �(φ 0 out ) b is the diffusivity outside droplets. Solving this equation and obtaining the corresponding flux j out outside the droplet is more difficult since the environment of the droplet might not be isotropic. We thus discuss the coupling of droplets to the dilute phase φ out in more detail in the next section.
If we know the fluxes j in and j out , we can determine the net accumulation of droplet material at the interface, which implies droplet growth. Note that only the normal components of the fluxes affect the shape, while the tangential components merely distribute material parallel to the interface. The shape changes of an isolated droplet are thus described by the interfacial speed v n in the normal direction 27 , see Supporting Information, Section I. General shape changes can result in non-spherical droplets, but since surface tension effects typically ensure a near-spherical shape, we project the general shape onto the degrees of freedom that we use to describe the droplet, where the integral is over the droplet surface, d is the space dimension, and S is the surface area of the droplet; see Supporting Information, Section II. Here, the first equation describes how material accumulates at the interface due to sum of all normal fluxes v n , leading to growth. The second equation describes how the weighted sum of all vectorial fluxes gives rise to drift. Taken together, Eq. (9) determines how an isolated droplet evolves in time. This involves Eqs. (4), (6), (8) as well as an approximation for the fluxes j out outside the droplet interface, which is the central part of our method that we discuss next.
Numerical model for many droplets. The dynamics of many droplets in the same system are coupled since they may exchange material via the dilute phase. To describe this exchange, and ultimately derive the flux j out at each droplet, we first consider the dynamics of the volume fraction of the dilute phase φ out . In principle, the dynamics of φ out follows from Eq. (7), with appropriate boundary conditions applied at the system's boundary and at all droplet surfaces. To simplify the description of the dilute phase, we assume that φ out is defined in the entire system, including where droplets are; see Fig. 1A. In this picture, droplets are local perturbations that exchange material with the background field φ out .
Dynamics of the background field. The background field φ out changes due to diffusion and reactions, even if droplets are absent. To capture this dynamics, we discretize the continuous field φ out on a uniform Cartesian grid with a distance x between the neighbouring support points. We then evolve Eq. (7) in time using finite central differences and explicit temporal stepping 45 . Since we do not need to resolve droplets at this scale, the spatial discretization can be much larger than in traditional Cahn-Hilliard equations.
Growth of a single droplet in a background field. To obtain the flux j out in the vicinity of an isolated droplet, we need to determine φ out in the region surrounding this droplet. We do this by considering Eq. (7) in an annular shell of thickness ℓ surrounding the droplet, which we further discretize into N sectors in the angular dimensions; see Fig. 1. In two-dimensional systems, we place sectors of equal size uniformly around the circle. Such symmetric placement is impossible in three dimensions, where we instead place N points approximately uniformly on the sphere and use a spherical Voronoi tessellation 46 to determine the corresponding sectors. In shell is the volume fraction of droplet material in the background field at the outer side of the m-th shell sector, which we estimate from a linear interpolation of the discretized background field φ out ; see Fig. 1. Since φ (m) out (r) typically varies only marginally in the shell sector, we also linearize the reac- . Taken together, we obtain an analytical approximation of φ (m) out (r) in each shell sector, from which we determine the local normal flux j (m) out outside the droplet; see Supporting Information, Section IV. Using these expressions together with Eqs. (6), (8) and (9), we find that individual droplets evolve according to where R and S are radius and surface area of the droplet, respectively, A m is the inner area of the shell sector and n m is the unit vector pointing from the droplet center to the m-th shell center; see Fig. 1. We use Eq. (10) to describe how internal reactions and external material exchange with the background affects the dynamics of each droplet.
Coupled dynamics of droplets and the background field. Equation (10) describe how droplets change when they exchange droplet material with the background field φ out . Due to material conservation, the material flux from the droplet toward each sector m, j (m) out · n A m , needs to accumulate in the background field. We use a linear interpolation at the midpoint of the inner boundary of the shell section (red points in Fig. 1B) to add the respective amount to the background field φ out . Note that negative fluxes j (m) out distribute material from the background field to the droplet and thus lead to growth. Taken together, this procedure ensures material conservation while preserving anisotropies of the exchange.
Full simulation. The full numerical method evolves the state of the system, i.e., the discretized background field φ out (r) and the positions x i and radii R i of all droplets, in time. We propose an explicit iteration, where the state at t + t is directly determined from the state at time t. Here, we first evolve the reaction-diffusion equation Eq. (7) of the background field and then iterate over all droplets. For each droplet, we determine the fluxes j (m) out · n for all shell sectors m and remove the associated material from the background field. We then update the droplet's position and radius according to Eq. (10). Starting from an initial state at t = 0 , this algorithm allows us to evolve the dynamics forward in time.
Choosing simulation parameters. The algorithm described above has several parameters that need to be chosen wisely for an accurate and fast simulation. In particular, we need to specify the discretization x of the background field, the shell thickness ℓ , the typical size s of a shell sector, and the time step t . We next discuss out at the outer side (black dot), which is determined from the background field using bilinear interpolation. The background field φ out is uniformly discretized on a Cartesian grid (gray grid).  (3). Our test case consists of two passive droplets of initial radius R 0 = 20w whose centers are separated by S d = 10R 0 . The droplets are placed in a background of vanishing initial volume fraction, φ out (r, 0) = 0 , so the system is under-saturated and the droplets will shrink. For the continuous model, we exploit the angular symmetry of the problem and consider an azimuthally symmetric cylindrical domain where r, z ∈ [0, 682w] . Conversely, our effective model is simulated in a 3-dimensional cubic domain of size [0, 1000w] 3 . We compare the final radii of the droplets after a duration T, where the droplets typically have shrunk by about 20% . The deviation of the mean droplet radius �R * � of our effective model compared to the radius R CM of the continuous model allows us to determine the crucial simulation parameters x , ℓ , and s.
Grid discretization x. The spatial discretization x determines the resolution at which variations of the background field φ out are resolved. Consequently, the choice of x is based on the problem: If spatial interactions are negligible and a mean-field model is desired, x can be arbitrarily large. Conversely, if spatial correlations between droplets are important, x needs to be smaller than the droplet separation. Another case are external gradients that affect droplets 12 , where x needs to be on the order of the droplet radii, so spatial anisotropies can be resolved on the droplet level.
Annular shell thickness ℓ. The most crucial part of our numerical method describes how material exchanges between the droplets and the background field. To describe this exchange faithfully, we interpolate the background field in an annular shell around the droplet. The thickness ℓ of this shell can thus be interpreted as an interpolation length scale and its value affects the accuracy of the simulation: If ℓ ≪ �x , the fluxes j (m) out are overestimated since they scale with ℓ −1 ; see Supporting Information, Section IV. Conversely, if ℓ ≫ �x , the background field would not be evaluated in the vicinity of the droplet, so interactions cannot be captured correctly. Taken together, we conclude that ℓ ∼ �x is a reasonable choice for the shell thickness. Indeed, Fig. 2A shows that this choice leads to a faithful estimate of the droplet growth for various values of x.
Shell sector width s. To resolve spatial anisotropies around a droplet, we discretize the shell into N sectors; see To obtain φ shell for each sector, the background field φ out is interpolated once per sector. Consequently, a larger number of sectors leads to a finer discretization and potentially a higher accuracy at the expense of larger computational cost. However, the accuracy is limited by the background discretization x as increasing N will have hardly any benefit if the distance s between interpolation points is already smaller than x . For two and three dimensions, we respectively use �s ≈ 2πR/N and �s ≈ 4πR 2 /N ; see Fig. 1B. Using s ∼ x , we can solve these equations for N, so that we can determine the number of sectors for each droplet based on its instantaneous radius R . Note that this implies that the dynamics of larger droplets will be described by more sectors to faithfully capture the interaction with their surrounding. Fig. 2B shows that the shell sector size s has only marginal effects in simple situations.   www.nature.com/scientificreports/ also render simulations unstable. We next separately analyze the dynamics of the background field, the shell, and the droplet growth to identify the maximal suitable value of t. The dynamics of the background are described by the partial differential equation Eq. (7), which we here solve using a simple explicit Euler scheme. A standard von Neumann stability analysis shows that this scheme is stable if �t < �x 2 /(2D out ) , where D out is the diffusivity in the background field. Consequently, a suitable time step for evolving the background field is �t background = 0.1�x 2 /D out , where we chose the constant pre-factor conservatively. Similarly, we define �t shell = 0.1ℓ 2 /D out for the shell. To ensure faithful dynamics of droplet growth, we demand that the relative growth R −1 |dR/dt| is small during a single time step t . Assuming that typical droplets are not much smaller than the mean initial droplet radius R , this implies a maximal time step �t drop = 0.1�R� 2 /D out . Finally, we also consider the time scale of reactions, �t reaction = 0.1/(max φ |s(φ)|) , based on the maximal rate of s(φ out ) . Taken together, we set the time step of the simulation to the minimal value of the four limiting time scales determined above.

Validation
We showed above that �x ≈ ℓ ≈ �s is a sensible choice for the parameters of our algorithm. To see how this choice affects accuracy and speed of the simulation, we next present three simulation scenarios, which range from single droplets in an heterogeneous environment to coarsening of large dilute emulsions.
Passive droplet in external gradient. We first consider a single droplet in an external composition gradient, which is maintained via boundary conditions. Biological cells use such a setup to control the position of droplets in their interior 10,12,47 . Fig. 3 shows that the effective droplet model captures the drift and growth of the passive droplet quantitatively. While the resulting dynamics are very similar, the run time of the simulations are very different: The continuous model took roughly one day to complete, while the effective model finished within 10 seconds on identical hardware. Since the continuous model is much slower, we performed some tests in the subsequent sections only with the effective model.
Active droplet with logistic growth. We next test whether our effective model also captures the growth of droplets subjected to non-linear chemical reactions. We here consider logistic growth, s(φ) = kφ(1 − aφ) , where k sets the reaction rate and a determines the chemical equilibrium; the parameter a −1 is often called the carrying capacity. In a phase separating system, these reactions produce droplet material outside the droplet (where φ < a −1 ) and destroy it inside (assuming φ 0 out > a −1 ). These reactions are thus qualitatively similar to the active droplets with linear chemical reactions that we discussed above. However, the non-linear chemical reactions now lead to a non-monotonous growth of the active droplet; see Mean-field coarsening of passive droplets. We next consider the interactions of many passive droplets in a dilute emulsion. When droplets only interact via the spatially averaged background field, Lifshitz and Slyozov predicted that the average droplet radius R grows as t 1/3 in this case 43,44,48 . Our simulation of 10 5 droplets indeed recovers this scaling (Fig. 5A) when we mimic this situation by setting the discretization x to the system size. Moreover, Fig. 5B shows that the distribution of radii also follows the universal shape Mean-field coarsening of active droplets. As a final example, we consider the interaction of many active droplets in a dilute emulsion. Here, we focus on the simplest case of a first-order reaction between the solvent and the droplet material, which is known to suppress Ostwald ripening 27,39 . We thus solve Eq. (3) using the reaction flux Fig. 6 shows that the emulsions with broad initial sizes quickly converge to mono-disperse distributions in 2 and 3 dimensions. The droplet size in stationary state is very close to the theoretical prediction, which we obtain numerically from the condition j in = j out using Eqs. (5) and (7). Taken together, we thus demonstrated that our effective model faithfully recovers important physical behaviour of active droplets.

Outlook and discussion
We showed that our effective method is orders of magnitudes faster than the continuous model while still accurately capturing the dynamics of droplets under the influence of chemical reactions and external gradients. To demonstrate that the method also extends to more challenging situations, we finally simulate the combination of chemical reactions and external gradients on the dynamics of the droplets. Fig. 7 shows that droplets grow as they drift along the gradient and they approach the fixed radius given by R 3D , so that this system controls droplet drift and size. Taken together, our simulations demonstrate that the novel simulation method captures the dynamics of interacting active droplets efficiently.
To gain the significant speed-up, our approach focuses on relevant degrees of freedom and leverages analytical results. This method can in principle be extended to more challenging situations in the future. For instance, droplets embedded in an elastic matrix affect each others growth, which can be described by similar effective theories 49,50 . Similar dynamics will also inform the dynamics of droplets in cells, where for instance chromatin  www.nature.com/scientificreports/ and the cytoskeleton suppress coalescence of condensates [51][52][53] . Our method can in principle also account for fluid flows, which are present in many liquid-like systems 10,54,55 . Here, large scale flows will advect the droplets and also affect the background field. Finally, we could account for Brownian motion of droplets, their coalescence upon contact, and their spontaneous division in sufficiently strongly driven systems 56,57 . In particular, droplets have shown anomalous coarsening behavior in cells, primarily due to hindrance and physical barriers which curb their ability of coalesce 51,52,58,59 . Extending our effective method to account for these physical processes will allow analyzing more and more complex situations, approaching the complexity necessary to understand the behavior of many droplets in biological cells.
In summary, we have demonstrated that our effective method faithfully captures the effects of chemical reactions and external chemical gradients. The method is several orders of magnitude faster than traditional continuous models, making it viable for fast and computationally efficient simulations of systems with many droplets. More importantly, our model provides a modular platform, which can be extended with other relevant physical phenomena affecting droplets, thus shedding insights on the formation, dissolution, stability, and sizes of biomolecular condensates.